clear
drop _all

cd "[FILE_LOCATION]" // Change to appropriate file location
set seed  123321


global sims 1000
global Ns 500 1000 2000
global Ds 0 1 2
global b 1
global cov .5,.1,.3,.1\.1,.5,.1,.3\ .3,.1,.5,.1\.1,.3,.1,.5

// Table E.1 - Continuous Instrument
mat Table = J(15,9,0)
capture program drop contcont
program define contcont, rclass
drop _all
set obs `1'

drawnorm eta0 eta1 ep0 ep1 , n(`1') cov($cov ) means(0,1,0,$b )
if `2'==0{
replace eta1 = eta0 + 1
}
if `2'==2{
	replace eta1 = eta0 + 1
    replace ep0 = .45*(exp(eta0)-exp((.5)/2)) + rnormal(0,.5)
	replace ep1 = $b + .2*(exp(eta0)-exp((.5)/2))+ rnormal(0,.65)
}
gen z = 2*runiform()
gen x = eta0 + eta1*z
gen y = ep0 + ep1*x

reg y x
return scalar b_ols = _b[x]
ivreg y (x = z)
return scalar b_iv = _b[x]
ivcrc y (x = z)
return scalar b_ivcrc = _b[x]
quietly reg x z
quietly predict vtilde, resid
quietly gen vtilde2 = (vtilde)^2
quietly reg vtilde2 z c.z#c.z
quietly predict vvar, xb
quietly replace vvar = 1e-6 if vvar<1e-6
quietly gen v = vtilde/vvar
reg y x v c.v#c.z c.v#c.x c.v#c.x#c.z
return scalar b_cf1 = _b[x]
reg y x vtilde c.vtilde#c.z c.vtilde#c.x c.vtilde#c.x#c.z
return scalar b_cf2 = _b[x]
end

foreach d of global Ds{
scalar scount = 1
foreach N of global Ns{
simulate b_ols = r(b_ols) b_iv = r(b_iv) b_ivcrc = r(b_ivcrc) b_cf1 = r(b_cf1) b_cf2 = r(b_cf2), reps($sims) : contcont `N' `d'
sum
quietly sum b_ols
mat Table[1,scount + 3*`d'] = r(mean) - $b
mat Table[2,scount + 3*`d'] = r(sd)
mat Table[3,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
quietly sum b_iv
mat Table[4,scount + 3*`d']= r(mean) - $b
mat Table[5,scount + 3*`d'] = r(sd)
mat Table[6,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
quietly sum b_ivcrc
mat Table[7,scount + 3*`d']= r(mean) - $b
mat Table[8,scount + 3*`d'] = r(sd)
mat Table[9,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
quietly sum b_cf1
mat Table[10,scount + 3*`d']= r(mean) - $b
mat Table[11,scount + 3*`d'] = r(sd)
mat Table[12,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
quietly sum b_cf2
mat Table[13,scount + 3*`d']= r(mean) - $b
mat Table[14,scount + 3*`d'] = r(sd)
mat Table[15,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
scalar scount = scount+1
}
}
mat list Table

putexcel set SimResult, modify sheet("Table E.1")
putexcel B1 = "Single FS Het."
putexcel E1 = "Double FS Het."
putexcel H1 = "Misspec. CLP"
putexcel A2 = "n="
putexcel B2 = "500"
putexcel C2 = "1000"
putexcel D2 = "2000"
putexcel E2 = "500"
putexcel F2 = "1000"
putexcel G2 = "2000"
putexcel H2 = "500"
putexcel I2 = "1000"
putexcel J2 = "2000"
putexcel A3 = "OLS"
putexcel A6 = "IV"
putexcel A9 = "IVCRC"
putexcel A12 = "CLP-CF"
putexcel A15 = "CLP-CF"
putexcel B3 = matrix(Table)



// Table E.2 - Binary Instrument
set seed  123321
global sims 1000
global Ns 500 1000 2000
global Ds 0 1 2
global b 1
mat Table = J(18,9,0)
global alpha00 = .45
global alpha01 = .95
global alpha10 = .75
global alpha11 = .45
global p = .6
global cov1 .5,-.25,.2,.1 \ -.25,.5,.1,.2 \ .2,.1,.5,.1\.1,.2,.1,.5
global cov2a .2,-.1,-.05,-.05 \-.1,.2,.1,.1 \ -.05,.1,.2,.05\-.05,.1,.05,.2
global cov2b .2,-.1,.05,.05 \-.1,.2,-.1,-.1 \ .05,-.1,.2,.05\.05,-.1,.05,.2
mat init = [0, ln((1-$p )/$p ), -.5, .75, -.3, .45,.2,.2]

capture program drop contbin
program define contbin, rclass
drop _all
set obs `1'

gen z = rnormal()
replace z = z>0

if `2'==0{
	drawnorm eta0 eta1 ep0 ep1 , n(`1') cov($cov1 ) means(-.5,1,0,$b )
	gen x = eta0 + eta1*z
	gen y = ep0 + ep1*x
}

if `2'==1{
	gen vgen = runiform()
	gen t0 = runiform()
	gen t1 = runiform()
	// z = 0, j = 0 
	gen a = -$alpha00 *(1 -2*vgen)
	gen b = (1+$alpha00 - 2*$alpha00 *vgen) 
	gen c =  - t0
	gen u00 = (-b + (b^2 - 4*a*c)^.5)/(2*a)
	drop a b c
	// z = 1, j = 0 
	gen a = -$alpha10 *(1 -2*vgen)
	gen b = (1+$alpha10 - 2*$alpha10 *vgen) 
	gen c =  - t0
	gen u10 = (-b + (b^2 - 4*a*c)^.5)/(2*a)
	drop a b c
	// z = 0, j = 1 
	gen a = -$alpha01 *(1 -2*vgen)
	gen b = (1+$alpha01 - 2*$alpha01 *vgen) 
	gen c =  - t1
	gen u01 = (-b + (b^2 - 4*a*c)^.5)/(2*a)
	drop a b c
	// z = 1, j = 1 
	gen a = -$alpha11 *(1 -2*vgen)
	gen b = (1+$alpha11 - 2*$alpha11 *vgen) 
	gen c =  - t1
	gen u11 = (-b + (b^2 - 4*a*c)^.5)/(2*a)
	drop a b c

	gen ep0 = (.5^.5)*((invgammap(5, z*u10 + (1-z)*u00)-5)/5^.5)
	gen ep1 = $b + (.5^.5)*((invgammap(5, z*u11 + (1-z)*u01)-5)/5^.5)
	gen vscaled = (.5 + .25*z)*invt(3, vgen)
	gen x = z + vscaled
	gen y = ep0 + ep1*x
}

if `2'==2{
	drawnorm a1 a2 a3 a4 , cov($cov2a ) means(-.3,-.5,-.3,-.5)
	drawnorm b1 b2 b3 b4 , cov($cov2b ) means(.45,.75,.45,.75)
	gen t = runiform()<$p
	gen eta0 = a1*t +b1*(1-t) - .5
	gen eta1 = a2*t +b2*(1-t) + 1
	gen ep0 = a3*t +b3*(1-t)
	gen ep1 = a4*t +b4*(1-t) + $b
	gen x = eta0 + eta1*z
	gen y = ep0 + ep1*x
}

reg y x
return scalar b_ols = _b[x]
ivreg y (x = z)
return scalar b_iv = _b[x]
ivcrc y (x = z)
return scalar b_ivcrc = _b[x]
quietly reg x z
// Simple CLP
quietly predict vtilde, resid
quietly sum vtilde if z == 1
gen vvar = (r(sd))^2
quietly sum vtilde if z == 0
replace vvar = (r(sd))^2 if z == 0
quietly gen v = vtilde/vvar
quietly gen vtilde2 = vtilde^2 - vvar
reg y x vtilde c.vtilde#c.z c.vtilde#c.x c.vtilde#c.z#c.x
predict y1, xb
return scalar b_cf1 = _b[x]
// FGM
if `2' == 1{
sort vtilde
quietly kdensity vtilde if z==1, at(vtilde) gen(v1pdf) nograph
quietly gen v1cdf = 0 if _n ==1 
quietly replace v1cdf = v1pdf*(vtilde - vtilde[_n-1]) + v1cdf[_n-1] if _n>1
quietly kdensity vtilde if z==0, at(vtilde) gen(v0pdf) nograph
quietly gen v0cdf = 0 if _n ==1 
quietly replace v0cdf= v0pdf*(vtilde - vtilde[_n-1]) + v0cdf[_n-1] if _n>1
quietly gen vcdf = v0cdf*(1-z) + v1cdf*z - 1/2
reg y x vcdf c.vcdf#c.z c.vcdf#c.x c.vcdf#c.z#c.x
return scalar b_cf2 = _b[x]	
}
else{
	return scalar b_cf2 =0
}
// Mixed Normal
if `2' == 2{
	capture noisily{
	fmm 2, emopts(iterate(1)): regress vtilde z
	mat init2 = e(b)
	mat init2[1,1] = init
	fmm 2, from(init2) emopts(iterate(100)): regress vtilde z
	}
if _rc==0{
gen mua = _b[vtilde:1.Class] +_b[vtilde:1.Class#c.z]*z
gen mub = _b[vtilde:2.Class] +_b[vtilde:2.Class#c.z]*z
scalar sda = _b[/var(e.vtilde)#1.Class]^.5 
scalar sdb = _b[/var(e.vtilde)#2.Class]^.5
estat lcprob, nose post
scalar p = _b[1.Class]
gen vpdf = p*normalden((vtilde - mua),sda ) + (1-p)*normalden((vtilde - mub),sdb )
gen V1 = normalden((vtilde - mua),sda )/vpdf
gen Vmix = (normalden((vtilde - mua),sda )- normalden((vtilde - mub),sdb ))/vpdf
gen V2 = normalden((vtilde - mub),sdb )/vpdf
gen V3 = V1*(vtilde-mua)/sda^2 
gen V4 = V2*(vtilde-mub)/sdb^2 
reg y x Vmix c.Vmix#c.x V3 c.V3#c.(z x) c.V3#c.z#c.x V4 c.V4#c.(z x) c.V4#c.z#c.x
return scalar b_cf3 = _b[x]	
}
	else{
		return scalar b_cf3 = .
	}
}
else{
		return scalar b_cf3 =0
}
end


global Ds 0 1 2
foreach d of global Ds{
scalar scount = 1
foreach N of global Ns{
simulate b_ols = r(b_ols) b_iv = r(b_iv) b_ivcrc = r(b_ivcrc) b_cf1 = r(b_cf1) b_cf2 = r(b_cf2) b_cf3 = r(b_cf3), reps($sims) : contbin `N' `d'
sum
quietly sum b_ols
mat Table[1,scount + 3*`d'] = r(mean) - $b
mat Table[2,scount + 3*`d'] = r(sd)
mat Table[3,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
quietly sum b_iv
mat Table[4,scount + 3*`d']= r(mean) - $b
mat Table[5,scount + 3*`d'] = r(sd)
mat Table[6,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
quietly sum b_ivcrc
mat Table[7,scount + 3*`d']= r(mean) - $b
mat Table[8,scount + 3*`d'] = r(sd)
mat Table[9,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
quietly sum b_cf1
mat Table[10,scount + 3*`d']= r(mean) - $b
mat Table[11,scount + 3*`d'] = r(sd)
mat Table[12,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5

if `d'==1{
quietly sum b_cf2
mat Table[13,scount + 3*`d']= r(mean) - $b
mat Table[14,scount + 3*`d'] = r(sd)
mat Table[15,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5

}
else if `d'==2{
quietly sum b_cf3
mat Table[16,scount + 3*`d']= r(mean) - $b
mat Table[17,scount + 3*`d'] = r(sd)
mat Table[18,scount + 3*`d'] = (r(sd)^2 + (r(mean) - $b )^2)^.5
}
scalar scount = scount+1

}
}
mat list Table

putexcel set SimResult, modify sheet("Table E.2")
putexcel B1 = "MVN"
putexcel E1 = "FGM"
putexcel H1 = "Mixed Normal"
putexcel A2 = "n="
putexcel B2 = "500"
putexcel C2 = "1000"
putexcel D2 = "2000"
putexcel E2 = "500"
putexcel F2 = "1000"
putexcel G2 = "2000"
putexcel H2 = "500"
putexcel I2 = "1000"
putexcel J2 = "2000"
putexcel A3 = "OLS"
putexcel A6 = "IV"
putexcel A9 = "IVCRC"
putexcel A12 = "CLP-CF"
putexcel A13 = "$V=\tilde{V}$"
putexcel A15 = "CLP-CF"
putexcel A16 = "$V$ in Eq. (\ref{eq: VFGM})"
putexcel A18 = "CLP-CF"
putexcel A19 = "$V$ in Eq. (\ref{eq: VFM})"

putexcel B3 = matrix(Table)
